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burden,  to  Department  of  Defense,  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports  (0704-0188),  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington,  VA  22202-4302. 
Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  any  penalty  for  failing  to  comply  with  a  collection  of  information  if  it  does  not  display  a  currently  valid 
OMB  control  number. 

PLEASE  DO  NOT  RETURN  YOUR  FORM  TO  THE  ABOVE  ADDRESS. 


1.  REPORT  DATE  (DD-MM-YYYY) 

September  2012 


2.  REPORT  TYPE 

Final 


3.  DATES  COVERED  (From  -  To) 


4.  TITLE  AND  SUBTITLE 

A  Description  of  the  Framework  of  the  Atmospheric  Boundary  Layer  Environment 
(ABLE)  Model 


5a.  CONTRACT  NUMBER 


5b.  GRANT  NUMBER 


5c.  PROGRAM  ELEMENT  NUMBER 


6.  AUTHOR(S) 

Yansen  Wang,  Chatt  Williamson,  and  Benjamin  MacCall 


5d.  PROJECT  NUMBER 


5e.  TASK  NUMBER 


5f.  WORK  UNIT  NUMBER 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

U.S.  Army  Research  Laboratory 
ATTN:  RDRL-CIE-D 
2800  Powder  Mill  Road 

Adelphi,  MD  20783-1197 _ 


8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 

ARL-TR-6177 


9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 


10.  SPONSOR/MONITOR’S  ACRONYM(S) 


11.  SPONSOR/MONITOR'S  REPORT 
NUMBER(S) 


12.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  unlimited. 


13.  SUPPLEMENTARY  NOTES 


14.  ABSTRACT 

This  technical  report  documents  the  framework  of  a  three-dimensional  Atmospheric  Boundary  Layer  Environment  (ABLE) 
model.  The  objective  of  the  project  is  to  develop  a  microscale  boundary  layer  meteorological  model  to  fill  the  gap  in  Army 
applications.  A  set  of  incompressible  equations  for  momentum,  energy,  and  scalar  is  chosen  specifically  for  the  microscale 
meteorological  modeling.  A  Cartesian  coordinate  with  a  structured  grid  is  selected  for  a  fast  computation  and  domain 
preprocessing.  A  finite  volume  numerical  method  with  a  collocated  grid  is  applied  in  this  framework.  Some  planned  features  for 
this  model  system,  such  as  turbulence  model  for  sub-grid  scale  and  the  immersed  boundary  method,  are  also  described  briefly. 
Bench  mark  test  cases  compared  with  the  laboratory  data  and  with  other  numerical  simulations  are  also  performed  and 
documented. 


15.  SUBJECT  TERMS 

Atmospheric  Boundary  Layer,  Finite  Volume  method.  Microscale  Meteorological  Model 


16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION 

OF 

ABSTRACT 

18.  NUMBER 

OF 

PAGES 

19a.  NAME  OF  RESPONSIBLE  PERSON 

Y ansen  W ang 

a.  REPORT 

Unclassified 

b.  ABSTRACT 

Unclassified 

C.  THIS  PAGE 

Unclassified 

uu 

32 

19b.  TELEPHONE  NUMBER  ( Include  area  code) 

(301)  394-1310 

Standard  Form  298  (Rev.  8/98) 
Prescribed  by  ANSI  Std.  Z39.18 


11 


Contents 


List  of  Figures  v 

1.  Introduction  1 

2.  Governing  Equations  2 

2.1  Mean  Atmospheric  V ariables . 2 

2.2  Turbulence  Closure  Models . 3 

2.2.1  k-s  model . 4 

2.2.2  Large  Eddy  S  imulation . 5 

3.  Numerical  Methods  for  the  Model  5 

3.1  The  Structured,  Collocated  Computational  Grid . 5 

3.2  Numerical  Method  for  Scalar  Transport  Equation . 7 

3.3  Numerical  Method  for  Momentum  Transport  Equations . 10 

3.3.1  Pressure- velocity  Coupling  in  the  Outer  Iteration . 1 1 

3.4  Method  for  Solving  the  Discretized  Algebraic  System  Equations . 14 

4.  Boundary  Conditions  15 

4.1  Inflow  Boundary . 15 

4.2  Outflow  Boundary . 15 

4.3  Wall  Boundary . 15 

4.3.1  Wind  and  Pressure . 15 

4.3.2  Temperature,  Thermal  Flux  from  Ground . 16 

4.4  Symmetry  Plane  (Free- slip) . 16 

4.5  Immersed  Boundary  Conditions . 17 

5.  Preliminary  Results  17 

5 . 1  Flow  Past  a  Infinite  Long  Cylinder . 17 

5 .2  Three-dimensional  Lid-driven  Cavity  Flow . 18 

6.  Conclusion  20 


iii 


7. 


References 


21 


List  of  Symbols,  Abbreviations,  and  Acronyms 


23 


Distribution  List 


24 


IV 


List  of  Figures 


Figure  1.  A  finite  volume  stencil  for  ABLE  model.  The  symbol  P  is  the  center  point  that  the 
FV  represents,  and  E,  W,  S,  N,  T,  B  represent  the  centers  of  the  neighboring  FVs.  The 
lowercase  letters  e,  w,  s,  n,  t,  b  represent  the  center  points  on  the  each  faces  of  the  FV  P. 

The  geometric  length  of  the  FV  (AX,  AY  AZ)  can  be  variable  for  non-uniform  structured 
computational  grid . 6 

Figure  2.  A  block  diagram  for  the  SIMPFE  algorithm . 14 

Figure  3.  A  schematic  diagram  for  the  wall  boundary . 16 

Figure  4.  A  schematic  diagram  for  top  symmetry  boundary . 16 

Figure  5.  A  comparison  of  the  ABFE  solution  with  the  laboratory  visualization  (Adopted 

from  Van  Dyke  1982)  at  different  Reynolds  numbers . 18 

Figure  6.  ABFE  model  simulation  of  a  lid-driven  cavity  flow.  The  cross  sections  denoted 
from  the  color-coded  planes.  The  results  were  compared  with  the  laboratory  test  data 
from  Prasad  and  Keosseff  (1989) . 19 

Figure  7.  ABFE  simulation  compared  with  other  numerical  model  results  for  the  YZ  plane  at 
X=0.75.  Feft  panel:  streamline  from  the  simulation  of  Zang  et  al.  (1994).  Right  panel: 
the  flow  field  cross  section  at  the  same  location  as  Zang  et  al . 20 


v 


Intentionally  Left  Blank. 


vi 


1.  Introduction 


In  meteorology,  microscale  phenomena  exhibit  spatial  scales  from  a  few  meters  to  hundreds  of 
meters  and  temporal  scales  from  minutes  to  several  tens  of  minutes.  Microscale  processes 
directly  impact  many  Army  applications — environmental  stress  factors  for  Soldiers’  operation; 
atmospheric  transport  and  diffusion;  electro-optical  (EO)  propagation;  and  operations  involving 
manned  and  unmanned  aircraft.  This  project  seeks  to  fill  an  operational  gap  in  weather 
forecasting  by  developing  a  numerical  model  uniquely  suited  for  forecasting  in  the  microscale, 
an  important  setting  for  the  Army.  This  report  details  the  initial  results  of  what  is  a  long  and 
significant  development  process. 

Currently,  there  is  no  microscale  community  meteorological  model  available.  With  the  objective 
of  accurate  and  timely  atmospheric  prediction  for  the  battlefield  and  military  installations,  we 
have  considered  some  alternatives  to  address  this  problem.  One  alternative  is  to  run  a  mesoscale 
model  with  microscale  temporal  and  spatial  resolution.  However,  state-of-the-art  numerical 
weather  prediction  (NWP)  models  are  designed  to  resolve  the  mesoscale  (from  one  to  tens  of  km 
and  one  to  several  hours)  up  to  regional  scale.  Mesocale  NWP  models  deal  with  larger  weather 
features  and  systems  such  as  fronts,  strong  storms,  and  precipitation.  Those  mesoscale  models 
cannot  simulate  fine-scale  wind  circulations  and  other  environmental  factors  due  to  their  neglect 
of  fine-scale  topography  such  as  small  hills,  buildings,  and  forests  in  the  atmospheric  boundary 
layer  (ABL).  The  physical  parameterizations  for  microscale  atmospheric  motion  are  also 
different  from  those  used  in  mesoscale  models  (Wyngaard,  2004;  Muschinsky  et  al.,  2004). 
Moreover,  as  NWP  models  move  to  finer  resolution  grids,  the  terrain-following  coordinate  in  the 
mesoscale  models  becomes  problematic.  The  vertical  coordinate  used  by  most  mesoscale  models 
is  not  designed  to  handle  the  sharp  variations  of  the  resolved  surface  features.  On  the  other  hand, 
currently  available  high-resolution,  computation  fluid  dynamics  (CFD)  models  are  mostly 
developed  for  specific  industrial  usages,  and  therefore,  lack  important  atmospheric  processes 
such  as  turbulent  transport  between  soil,  urban,  vegetation,  and/or  surface  water  and  the 
atmosphere  and  radiation.  We  seek  to  fill  a  gap  in  Army  capabilities  by  developing  a  prognostic, 
physics-based,  microscale  model  optimized  for  use  in  the  ABL.  Furthermore,  we  hope  to 
transition  the  Atmospheric  Boundary  Layer  Environment  (ABLE)  model  to  a  real-time 
prediction  tool  as  future  computers  become  more  powerful. 

In  last  several  years,  we  have  developed  a  three-dimensional  wind  field  model  (3DWF)  for 
application  over  complex  terrain  and  in  urban  settings  (Wang  et  al.  2005;  Wang  et  al.  2010; 
Hanna  et  al.  2011).  Although  a  diagnostic  type  model  such  as  3DWF  is  reasonably  accurate  for 
strong  and  neutral  wind  conditions  and  consumes  very  little  computing  power,  it  does  not  model 
other  needed  atmospheric  variables  (i.e.,  temperature,  moisture,  and  radiation)  and  has  intrinsic 
shortcomings  in  dealing  with  convective  and  stable  conditions.  For  these  reasons,  we  have 
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initiated  a  new  project  in  FY2012  to  develop  a  next  generation  microscale  meteorological  model, 
the  ABLE  model,  to  meet  the  Army  application  requirements  for  characterizing  the  microscale 
motion  of  the  atmospheric  boundary  layer.  The  ABLE  model  is  planned  to  be  an  advanced 
microscale  meteorological  model  that  couples  the  multiple  atmospheric  environmental  variables 
such  as  wind,  temperature,  moisture,  and  scalars  or  pollutant  transports.  The  model  is  based  on  a 
set  of  three-dimensional,  prognostic,  incompressible,  Navier-Stocks  equations  with  special 
attention  given  to  the  microscale  variation  of  atmospheric  variables  caused  by  complex 
boundaries  such  as  buildings,  forest  canopies,  and  complex  terrain.  The  model  will  take 
advantage  of  recent  developments  in  the  treatment  of  boundary  conditions,  computational  grid 
generation,  numerical  solvers,  and  turbulence  modeling.  This  model  will  be  a  large  leap  forward 
in  our  capability  to  model  the  microscale  ABL  variables  and  their  interactions. 


2.  Governing  Equations 


2.1  Mean  Atmospheric  Variables 

The  proposed  new  microscale  model  will  focus  on  the  flows  in  the  ABL.  There  are  several 
special  considerations  for  a  microscale  meteorological  model.  The  framework  for  the  ABLE 
model  is  based  on  the  following  conservation  equations  for  the  meteorological  variables, 
including  wind,  temperature,  pressure,  scalars  such  as  moisture  or  pollutant  concentrations,  and 
energy  transfer  between  the  Earth’s  surface  and  the  atmosphere.  We  use  a  set  of  incompressible 
Navier-Stocks  system  equations  with  the  Boussinesq  approximation  to  eliminate  sound  waves 
that  have  no  meteorological  significance.  The  advantage  of  not  resolving  sound  waves  is  that  the 
model  time  marching  scheme  will  allow  a  much  larger  time  step  according  to  the  Courant- 
Friedrichs-Lewy  (CFL)  criteria.  The  Coriolis  force  is  rationally  neglected  because  the  domain 
that  this  model  intends  to  cover  is  very  small  and  the  Coriolis  force  is  several  orders  of 
magnitude  smaller  than  other  forces.  The  choice  of  this  equation  set  is  justified  because  we  are 
modeling  very  fine-scale  atmospheric  motions  in  the  ABL  where  deep  moist  convection  is  not 
present  (Stull  1989;  Durran  2008).  The  first  stage  in  the  development  is  focused  on  the  Reynolds 
Averaged  Navier-Stocks  (RANS)  type  model,  and  the  second  stage  of  the  development  is 
extended  to  the  Large  Eddy  Simulation  (LES)  type  of  model.  In  derivation  of  the  following 
model  equations,  two  operations  were  performed  for  the  meteorological  variables.  The  first  step 
is  to  represent  the  variables  of  pressure,  temperature,  and  density  as  the  sum  of  their  base  states 
and  deviations.  The  density  deviation  from  inertial  terms  are  all  neglected  except  in  the  gravity 
terms,  and  the  ratio  of  density  deviation  to  the  base  density  is  then  replaced  with  the  ratio  of 
potential  temperature  deviation  to  the  base  state  potential  temperature  in  the  gravity  term  (i.e., 
Boussinesq  approximation).  The  second  step  separates  the  velocity,  potential  temperature,  and 
the  scalar  variables  into  the  two  parts,  the  mean  and  the  turbulent,  by  ensemble  averaging  (for 
RANS)  or  spatial  filtering  (for  LES).  The  detailed  derivations  of  the  model  equation  set  can  be 
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found  in  reference  books  (Stull  1989;  Wyngaard  2010).  The  following  is  the  RANS  equation  set 
that  is  written  in  conservative  tensor  form  in  the  Cartesian  coordinate: 


d(pU,)  dipUJJj)  dP  d2Ui  g0  dipu'u)) 

- —  H - —  = - h  V - +  On - — 

dt  dxj  dxi  dx2j  90  dxj 

d(p9)  djpUjd )  _  d20  1  QRn  EE  dipu'jO') 

dr  dxj  9  dx2  Cp  dxj  Cp  dxj 

d(pC)  d(pU  jC)  d2(pC)  dipu^) 

- 1 - - -  =  Vr  — r - - - h  S 

dt  dxj  dx2  dxj 

dT 

Cg-^  =  Rn-LvE-H  +  Hm 


(1) 

(2) 

(3) 

(4) 

(5) 


where  I/,  and  U  .  are  the  ensemble  averages  of  wind  velocity  components;  prime  terms  are  the 

turbulent  fluctuations  of  the  variables;  over-bar  terms  are  the  turbulent  flux  or  Reynolds  stress 
terms;  p  and  90  are  the  base  references  of  air  density  and  potential  temperature;  C  is  the  mean 

scalar  mass  mixing  ratio  and  the  scalar  that  can  be  moisture  or  other  gases;  and  P  and  9  are  the 
pressure  and  potential  temperature  deviations  from  the  base  state.  The  averaged  turbulent  flux 
terms  of  momentum,  heat,  and  scalar  parameters  are  represented  by  the  over  bar  of  primed  terms. 
Tg,  R„,  E,  H,  Hm  are  the  ground  surface  temperature,  net  radiation  flux,  evaporation,  sensible 
heat,  and  molecular  heat  conduction  between  the  atmosphere  and  the  ground  surface, 
respectively;  v,v0,vcare  the  dynamic  molecular  viscosity,  diffusivity  of  heat,  and  diffusivity  of 

scalar  parameters,  respectively;  Cp,Cg,Lv,S  are  the  specific  heat  of  air  at  constant  pressure,  the 

heat  capacity  of  the  Earth  surface  materials,  the  latent  heat  of  evaporation,  and  source  strength 
term  for  scalar  parameters,  respectively;  and  Sj3  is  the  Kronecker  delta  for  tensor  notation.  The 

conservation  of  mass,  momentum,  energy,  scalar  parameters,  and  radiation  flux  are  governed  by 
the  equations  1-5,  respectively,  in  tensor  differential  form.  Note  that  absolute  pressure  is  of  no 
significance,  only  the  pressure  gradient  affects  the  flow  in  the  incompressible  flow.  The  pressure 
field  is  governed  by  a  Poisson  equation,  which  can  be  derived  by  using  the  pressure-velocity 
coupling  to  satisfy  the  continuity  equation. 

2.2  Turbulence  Closure  Models 

There  are  many  turbulence  closure  models  that  have  been  developed  for  different  flows.  In  this 
section,  we  only  introduce  the  standard  k-  s  model  and  the  LES  model.  This  subject  will  be 
described  in  detail  once  we  complete  the  implementation  of  a  turbulent  parameterization  in  the 
ABLE  model. 


3 


2.2.1  k-s  model 

In  the  RANS-type  turbulence  model,  the  turbulent  fluxes  result  from  partitioning  total  velocity 
variables  into  fluctuation  and  mean  parts  via  either  time  or  ensemble  averages.  The  turbulent 
fluxes  are  computed  using  the  turbulent  viscosity,  which  is  analogous  to  molecular  viscosity  in 
the  laminar  flow.  In  the  k-s  model  (Launder  and  Spalding  1974),  the  turbulent  viscosity  (Km)  is 
determined  by  the  turbulent  kinetic  energy  (k  =  u’u’  12)  and  dissipation  rate  ( s )  by 


K 


m 


(6) 


where  Cu  is  an  empirical  coefficient.  The  turbulent  viscosities  for  heat  and  scalar,  Kh ,  and  Kc , 
are  related  to  the  K  via  the  effective  Prandtl  number  (Pr  =  K  /  K,  =  0.7)  and  Schmidt  number 
(Sc  =  Km  / Kc  =  0.9).  The  turbulent  kinetic  energy  (k)  and  turbulent  dissipation  rate  (s)  are 
computed  using  the  following  prognostic  equations: 
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(8) 


where  the  numerical  values  of  closure  coefficients  are  taken  as  follows: 

Csl=  1-44,  Ce2  =  1.92,  CM  =0.09,  ak  =1.0,  ^  =1.3 

Launder  and  Spalding  (1974)  were  the  first  ones  that  constructed  a  numerical  model  and  used 
above  closure  coefficients.  Detail  derivation  of  above  k  and  s  equations  can  be  found  in  Wilcox 
(2006).  The  physical  interpretations  of  the  terms  in  the  turbulent  kinetic  energy  (TKE)  equation  7 
are  as  follows.  The  terms  on  the  left  hand  are  the  tendency  and  advection  of  TKE  by  mean  flow. 
The  first  and  last  terms  on  the  right-hand  side  represent  the  shear  and  buoyancy  production  of  the 
TKE,  and  the  third  term  represents  the  transport  of  TKE  due  to  molecular  diffusion  and 
turbulence. 


The  turbulent  fluxes  terms  in  the  governing  equations  2-4  are  parameterized  with  the  resolvable 
mean  variables  using  the  Boussinesq’s  turbulent  viscosity  approach: 


U'U'j  =  2Sijk~K» 
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U':C  =  ~K 

1  dxj 

where  k  is  the  turbulent  kinetic  energy.  Km ,  Kh ,  and  Kc  are  the  turbulent  (eddy)  viscosities  of 
momentum,  heat,  and  scalar,  respectively. 

2.2.2  Large  Eddy  Simulation 

In  a  LES  model,  the  variables  are  partitioned  into  a  grid  resolvable  part  and  a  spatial  filtered 
fluctuation  part.  The  LES  is  a  modeling  technique  to  resolve  the  large,  energy  containing 
turbulent  eddies  and  the  subgrid  turbulent  fluxes  are  parameterized  using  either  an  eddy  viscosity 
model  or  a  dynamical  model  (Smagorinsky  1963;  Germano  et  al.  1990).  The  model  grids  have  to 
be  smaller  than  the  large  eddies  to  resolve  them.  Since  the  scales  of  the  large  eddies  in  many 
ABL  flows  are  quite  small,  the  computation  grids  are  much  finer  than  RANS-type  modeling.  The 
LES  requires  much  more  computational  power  compared  with  the  RANS.  We  defer  description 
of  the  LES  modeling  to  the  future  in  the  ABLE  model  development  processes.  Detail  discussion 
about  LES  can  be  found  in  Pope  (2000)  and  Wyngaard  (2010). 


3.  Numerical  Methods  for  the  Model 


3.1  The  Structured,  Collocated  Computational  Grid 

A  box  (hexahedron)  shaped,  structured  computational  grid  is  employed  in  the  ABLE  model.  The 
rationale  to  choose  the  structured  box-type  grid  is  that  generation  of  this  type  of  grid  is  much 
simpler  than  generating  the  unstructured  body-fitted  grid.  With  recent  research  and  development 
of  immersed  boundary  (IB)  methods  (Mittal  and  Iaccarino  2005;  Tseng  and  Ferziger  2003),  the 
complex  boundary  can  be  treated  with  the  box-type  structured  grid,  which  significantly  reduces 
computation  time  and  complexity.  The  body-fitted  computational  grid  generation  is  a  laborious 
process  often  taking  several  days  for  a  complex  terrain  (e.g.,  small  urban  setting),  but  the  box- 
type  structured  grid  would  generated  within  a  minute  for  the  same  computational  domain.  There 
are  other  incentives  for  using  the  box-type  structured  grid  in  the  computation,  including  (1)  less 
computational  time  required  for  one  iteration  on  the  structured  grid  in  comparison  with  the 
unstructured  grid;  and  (2)  the  stability  of  the  numerical  algorithm  on  the  structured  grid  being 
essentially  higher  than  on  the  unstructured  grid  (Tseng  and  Ferziger  2003). 

A  finite  volume  (FV)  method  is  used  to  discretize  the  model  equations.  The  FV  has  an  advantage 
over  finite  difference  methods  for  satisfying  mass  conservation  automatically  without  special 
treatments.  In  this  model,  the  FV  sizes  are  variable  to  accommodate  the  resolution  requirement  in 
different  situations.  For  an  example,  much  denser  grids  might  be  needed  to  deal  with  the  large 
gradient  of  environmental  variables  to  achieve  the  desired  accuracy  near  the  ground  surfaces  or 
buildings.  Figure  1  shows  an  arbitrary  control  volume  used  to  derive  the  discretized  numerical 
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equations  for  the  model.  The  center  of  this  FV  is  denoted  by  the  point  P,  and  the  uppercase 
letters,  E,W,N,S,T,  and  B  represent  the  neighboring  FVs’  centers  at  east,  west,  north,  south,  top, 
and  bottom,  respectively.  The  lowercase  letters,  e,w,n,s,t,  and  b  represent  the  locations  that  are 
half-distance  points  from  the  P  to  neighboring  FVs’  centers.  The  ABLE  model  also  uses  a 
collocated  grid  for  different  environmental  variables,  i.e.,  all  the  variables  are  solved  in  the  same 
location  denoted  as  P  for  a  control  volume.  The  collocated  grid  offers  various  advantages 
compared  with  the  traditional  staggered  grid  arrangement.  This  collocated  grid  setup  saves 
computational  time  and  memory  by  using  the  same  geometric  parameters  for  the  different 
variables.  The  collocated  grid  allows  for  a  single  set  of  common  grid  geometric  parameters 
reducing  computational  and  memory  requirements.  It  also  offers  great  deal  of  convenience  in 
boundary  treatment  for  different  variables  since  all  variables  coincide  with  the  boundary  of 
physical  boundaries  of  the  solution  domain.  The  property  of  the  collocated  grid  offers  a  great 
deal  of  simplicity  in  treatment  of  boundary  when  the  IB  method  is  used.  The  collocated  grid 
requires  a  special  interpolation  scheme  (Rhie  and  Chow  1983)  to  prevent  the  oscillating  velocity 
and  pressure  fields,  which  is  described  later. 


Figure  1.  A  finite  volume  stencil  for  ABLE  model.  The  symbol  P  is  the 
center  point  that  the  FV  represents,  and  E,  W,  S,  N,  T,  B 
represent  the  centers  of  the  neighboring  FVs.  The  lowercase 
letters  e,  w,  s,  n,  t,  b  represent  the  center  points  on  the  each  faces 
of  the  FV  P.  The  geometric  length  of  the  FV  (AX,  AY  AZ) 
can  be  variable  for  non-uniform  structured  computational  grid. 
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3.2  Numerical  Method  for  Scalar  Transport  Equation 


The  scalar  transport  equation  is  used  to  demonstrate  the  numerical  discretization  without  losing 
the  generality  for  other  governing  equations,  since  the  momentum  and  heat  transport  equations 
are  very  similar  to  scalar  transport  except  the  pressure  and  source  terms.  The  treatment  of  the 
pressure  gradient  force  in  the  momentum  equation  is  described  separately.  Following  Patankar 
(1980),  the  scalar  transport  equation,  including  the  turbulence  transport  parameterization,  can  be 
written  in  following  form: 


d(pC) 
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where  U,  V,  W  are  three  Cartesian  components  of  the  wind  velocity  and  S  is  the  source  strength 
term  for  the  scalar.  Using  an  implicit  time  differencing  scheme,  the  integration  of  this  equation 
over  a  box-shaped  control  volume  (figure  1)  gives 
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where  superscript  m  denote  the  values  from  last  outer  iteration  (time  marching),  Ax,  Ay,  Az  arc 

the  finite  control  volume  side  lengths  in  the  x,  y,  and  z  directions,  which  are  prescribed  at 
different  spatial  locations  for  a  non-uniform  grid  configuration.  Mass  fluxes  0cv  /acas  through 
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These  terms  incorporate  both  the  advection  and  diffusive  fluxes  through  the  each  of  the  control 
volume’s  faces.  Note  that  the  mass  flux  pU  uses  the  U  values  from  last  outer  iteration  m.  The 
mass  mixing  ratio  of  scalar,  C,  is  at  the  time  step  m+1  since  the  implicit  scheme  is  applied  to 
allow  the  large  time  steps  without  numerical  instability  problem.  The  superscript  of  the  time  step 
m+1  is  abbreviated  for  clarity  in  the  derivation.  The  continuity  equation  is  satisfied  in  the  scalar 
transport  equation.  The  advection  air  mass  fluxes,  F,,  are  determined  by  the  upstream  mass 
fluxes  (i.e.,  Fe  =  pU"'AyAz )  using  the  velocity  values  from  last  outer  iteration,  named  as 

upstream  difference  scheme  (UDS).  The  upwind  mass  flux  for  scalar  is  denoted  in  QUDS: 

QTS  =  p(UC)eAyAz  =  CP  max(  Fe, 0)  -  CE  max( -F„0)  =  CP  max(  Fe, 0)  +  CE  min(F„0) 

QTS  =  p(UQwAyAz  =  Cw  max(  Fj 0)  -  Cp  max( -F  ,0)  =  Cw  max(  F  ,0)  +  CP  min(F  ,0) 

QT  =  p(VC)n  AxAz  =  CP  maxi  F,  ,0)  -  CN  max(  -F„  ,0)  =  CP  max(  F,  ,0)  +  CN  min(  Fn  ,0) 

Q':ns  =  p(VC)s  AxAz  =  Cs  max(  F  ,0)  -  CP  max(  -F  ,0)  =  C,  max(  F  ,0)  +  CP  min(  F  ,0) 

(2,um  =  p(WC),  ArAy  =  C,,  max(  F,  ,0)  -  CT  max(  -Ff  ,0)  =  CP  max(  F,  ,0)  +  CT  min(  F,  ,0) 

QTS  =  p(WC)hAxAy  =  Cr:  max(  F„0)  -  C,,  max(  -F„0)  =  C/;  max(  Ffc,0)  +  CP  min(F„0)  (15) 

To  achieve  second-order  accuracy,  we  use  the  deferred  correction  approach  (Khosla  and  Rubin, 
1974,  Ferziger  and  Peric,  2002)  to  compute  the  advective  mass  flux  for  scalar,  which  includes 
the  central  difference  scheme  (CDS).  For  the  mass  flux  through  the  face  e  of  CV, 

Qe=Qr+(QeDS~QeDS)m-  (16) 

The  superscript  m  again  means  that  computation  of  those  fluxes  using  the  last  outer  iteration 
solved  value.  The  outcome  of  the  outer  iteration  is  that  UDS  contribution  is  canceled  out,  leaving 
a  CDS  solution  which  is  the  second-order  accuracy.  The  CDS  computation  of  scalar  mass  flux  is 
defined  as 


QeDS  =  QEre  +  QPd-re )  (17) 

where  r  =  (xe  -xp)/(xE  —xp) ,  QE  =  pUECEAyAz  ,  and  the  lowercase  subscripts  represents  the 
advective  fluxes  at  the  corresponding  faces  of  control  volume  and  the  uppercase  subscript  of  C 
denote  the  concentration  values  at  the  control  volume  centers.  The  concentration  gradients  in 
each  direction  are  discretized  in  following  central  differentials  for  the  computation  of  diffusive 
fluxes: 
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Here  we  only  show  cDe  -Oiv  flux  in  x  direction,  the  fluxes  in  other  two  (y  and  z.)  directions  are 
similar.  Let  D  =  p(yc  +  Kc)  be  the  total  diffusion  coefficient  for  a  scalar,  then 
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(19) 


Adding  the  diffusive  and  advective  fluxes  in  y  and  z  directions  and  after  some  simple  algebra 
manipulations,  the  discretized  scalar  transport  equation  can  be  written  as  following  liner  system 
equations: 

dpCp  +  dp-Cp  +  dwCw  +  dNCN  +  d^Cg  +  djCp  +  dgCjj  =  b  +  Q ,  (20) 


where 
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aB  =  -max(Ffc,0) 


(yc  +  Kc) 

Zp  ~  Z-B 


AxAy 


b  =  SAxAyAz  +  a‘P  C'P 


ap  —  ap  aE  aw  aN  as  aT  aB+(Fe  Fw )  +  ( Fn  Fs~)  +  (Ft  Fh ) 

At  this  point,  the  continuity  equation  1  can  be  written  in  similar  finite  volume  form: 

Fe-Fw  +  FN-Fs+FT-FB=0  (21) 

Therefore,  combining  the  equations  20  and  21,  we  have 

dp  dyy  d  yj  dp  flg  (22) 

In  equation  20,  Q=Qe+Qw+Qn+Qs+Qt+Qb  is  the  sum  of  source  terms  of  advection  in  each 
direction  due  to  the  deferred  correction  when  using  old  wind  speed  values  from  the  last  iteration 
(t).  (The  diffusion  source  term  is  neglected.)  For  example,  the  deferred  flux  at  e  of  CV  is 
expressed  as 

Qe  =  KQe)UDS -(Qe)CDST  (23) 

Deferred  mass  fluxes  at  other  faces  of  the  CV  can  be  computed  in  the  same  way. 

3.3  Numerical  Method  for  Momentum  Transport  Equations 

Comparing  equation  2  with  the  scalar  transport  equation  4,  besides  the  terms  of  change  of 
momentum  with  time,  advection  and  diffusion  in  both  equations,  there  are  two  more  terms — the 
pressure  gradient  and  the  gravity  in  the  momentum  equations.  The  momentum  equations  2  differ 
from  passive  scalar  transport  equation  4  in  two  important  ways:  the  momentum  equation  is 
nonlinear  in  advection  terms  and  coupled  between  all  the  velocity  components  so  that  the 
equation  must  be  solved  iteratively;  and  the  pressure  gradient  forces  require  solving  the  pressure 
field.  The  discretization  for  the  advection  and  the  diffusion  terms  are  the  same  as  for  the  scalar 
transport.  The  pressure  gradient  and  gravity  terms  need  special  treatment.  The  discretized 
algebraic  equation  for  equation  2  is  similar  to  the  scalar  transport  equation  except  the  extra 
pressure  and  buoyancy  terms.  The  velocity  components  from  last  outer  iteration  are  also  used  to 
determine  the  convective  and  diffusive  mass  fluxes.  Without  repeating  much  of  the  derivation  for 
the  scalar  equation,  the  discretized  momentum  equations  can  be  written  as  follows  for  U,  V,  and 
W  components: 

op* 
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«VX  =  IX,,  C -^-AxAyAz  +  SV 
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r)P 

o>;  =  W,\  -  —  AxAyAz  +  SW 

OZ 


(24) 


Note  that  in  above  algebraic  equation,  we  used  superscript  for  U,  V,  W.  and  P.  This  notation 
is  used  for  the  convenience  of  the  description  of  the  pressure- velocity  coupling  in  the  iteration 
algorithm  discussed  in  the  next  section.  The  nb  represents  all  the  six  neighbor  points,  E,  W,  N,  S, 
T.  B  (figure  1)  surrounding  to  the  point  P.  The  SU,  SV,  SW  are  source  terms.  Since  we  use  the 
collocated  grids,  the  pressure  term  need  to  be  interpolated  to  the  location  of  the  control  volume 
faces  as  denoted  in  the  above  equation.  Recall  that  only  the  deviations  of  the  pressure  from  the 
base  state  contribute  to  the  momentum  transport,  as  stated  in  basic  equation  section.  The 
computation  for  those  coefficients  is  the  same  as  in  the  discretized  scalar  equation.  The  last  terms 
on  the  above  equations  are  the  source  terms  including  the  buoyancy,  turbulent  diffusion,  and 
deferred  correction  in  the  inner  iteration. 


3.3.1  Pressure-velocity  Coupling  in  the  Outer  Iteration 

In  the  incompressible  equation  set,  there  is  no  independent  equation  for  pressure 
(thermodynamic  pressure  is  undefined).  Instead,  the  pressure  gradient  forcing  is  calculated  by 
taking  the  divergence  of  the  momentum  equation  2,  using  the  incompressible  continuity  equation 
1,  and  assuming  the  constant  density  and  viscosity  to  produce  the  following  Poisson  equation  for 
pressure: 

d2P  _  d 
dxf  dxi 


d(pUJUi) 


dXj 


(25) 


Note  that  the  absolute  pressure  value  is  not  needed;  only  the  gradient  of  the  pressure  is  important 
in  incompressible  flow.  The  pressure  deviation  contributes  to  the  pressure  gradient  force.  The 
pressure  and  velocity  are  not  solved  simultaneously;  rather,  they  are  solved  by  an  iterative 
process.  In  the  Semi-Implicit  Method  for  Pressure-Linked  Equations  (SIMPLE)  iteration 
(Patankar,  1980),  pressure  and  velocity  are  corrected  by  adding  the  correction  parts  (denoted  by 
superscript ')  for  the  variables  solved  in  last  outer  iteration  (denoted  by  superscript  *),  i.e., 


P  =  P*  +  P' 

U  =  U*  +  U' 
V  =V*+V' 
W  =W*  +W' 


(26) 


The  corrected  velocities  can  be  similarly  expressed  as  equation  24: 
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dz 


(27) 


Subtracting  the  equations  24  from  equations  27,  an  equation  set  for  correction  terms  is  derived: 


a 


r)P' 

lU'p  =YJ°nb  U'nb  ~~^AxAyAZ 


aVPV'p  =  2X*  Kb  -  —  AxAyAz 

oy 

OpW'p  =  YjCIZ,  Kb  -  —^—AxAyAz 

dz 


(28) 


The  SIMPLE  procedure  estimates  the  above  equations  by  omitting  K/K,  U  'nb  >  Kz, .  and 

Kb terms-  Since  this  is  an  iteration  procedure,  it  is  a  valid  approximation.  The  detailed 

discussion  can  be  found  in  Patankar  (1980)  and  Tu  et  al.  (2008).  Substituting  the  equations  28 
without  those  X  terms,  the  simplified  equation  can  be  written  as 
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where  B  =  AxAyAz, / aP ,  B  =  AxAyAzI aP,  B  =  AxAyAzI aP  .  By  summing  up  the  three 
equations  in  equation  29,  taking  the  divergence  of  the  resulted  summation,  and  using  the 
continuity  equation,  the  equation  for  P'  is  derived: 
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(30) 


This  Poisson  equation  can  be  discretized  in  similar  way  as  done  for  the  scalar  and  momentum 
equation. 

One  thing  that  needs  to  be  stressed  here  is  the  treatment  of  pressure  and  velocity  coupling  on  the 
collocated  grid.  The  Rhine  and  Chow  (1983)  velocity  interpolation  is  used  for  the  wind 
components  at  the  CV  faces  to  eliminate  spurious  pressure  and  velocity  oscillations  from  the 
collocated  grid.  In  contrast  to  the  simple  interpolation  of  velocity  components  only,  the  Rhine 
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Chow  velocity  interpolates  both  the  velocity  components  and  the  pressure  gradient  for  their  CV 
face  values,  and  adds  a  higher  order  correction  term  for  CV  face  velocity  interpolation.  This  is 
equivalent  to  adding  a  pressure  smoothing  term  that  is  proportional  to  the  third  derivative  of  the 
pressure  to  prevent  from  occurrence  of  spurious  pressure  mode.  This  interpolation  can  be 
expressed  as  follows  for  a  one-dimensional  example: 


1  SV 
Ue=~(UE+UP)  +  — 

2  2  aD 
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dx 
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dP_ 

dx 
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(31) 


The  SIMPLE  algorithm  is  summarized  in  figure  2.  The  initial  step  begins  by  guessing  a  pressure 
P*,  the  velocity  components,  U\  V*,  W  can  be  solved  using  equation  24.  The  pressure  (i3’)  field 
then  can  be  solved  via  Poisson  equation,  equation  30.  At  the  next  step,  the  pressure  and  velocity 
are  corrected  using  equation  26.  At  this  point,  by  using  the  corrected  velocity  components  and 
pressure,  the  temperature,  other  scalar  variables,  and  turbulence  transport  equations  can  be 
solved  using  the  velocity  and  pressure  values  from  the  last  step.  The  last  step  is  to  check  the 
convergence.  If  the  convergence  is  not  achieved,  replace  the  U* ,  V ,  W*  P*  with  the  newly 
solved  U,  V,  W,  P,  and  go  through  another  outer  iteration. 
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Figure  2.  A  block  diagram  for  the  SIMPLE  algorithm. 

3.4  Method  for  Solving  the  Discretized  Algebraic  System  Equations 

After  FV  discretization,  each  of  the  governing  equations  of  wind  components,  energy,  and 
scalars  becomes  a  set  of  linear  system  equations,  whose  coefficient  is  a  sparse,  seven-diagonal 
matrix.  There  are  several  methods  to  solve  this  type  of  system  of  algebraic  equations,  the 
successive  over  relaxation  (SOR)  method;  the  alternate-direction  implicit  (ADI)  method;  the 
strongly  implicit  procedure  (SIP)  (Stone,  1968;  Leister  and  Peric,  1994);  and  the  bi-conjugate 
gradient  stabilized  (BI-CGSTAB)  (Van  den  Vorst,  1992).  Since  most  computation  is  in  this 
matrix  solver,  we  will  continually  evaluate  solution  methods  for  the  ABLE  model  based  on 
speed,  accuracy,  and  numerical  stability;  some  experimentation  will  be  required  during  model 
development.  Terrain,  building,  and  morphology  conditions  may  also  alter  the  performance  of 
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the  solvers.  For  example,  a  technique  might  perform  better  in  urban  but  not  mountainous 
complex  terrain.  We  intend  to  perform  sensitivity  tests  for  these  numerical  solvers  and  provide 
guidance  to  the  user  for  which  options  to  use  when  setting  up  a  model  run.  For  the  prototype,  we 
choose  the  SIP  method  for  testing  the  framework.  The  SIP  method  creates  an  easily  factorable 
matrix  operator  for  iteration.  The  rate  of  convergence  is  increased  by  adding  a  contribution  based 
on  a  smoothness  assumption  in  the  dependent  variable. 


4.  Boundary  Conditions 


The  boundary  conditions  discussed  in  this  section  are  applied  only  for  the  mean  conservation 
equations  of  momentum,  energy,  and  scalar.  The  discussion  of  boundary  conditions  for  turbulent 
fluxes  and  TKE  equation  will  be  postponed  until  future  development.  We  also  planned  to 
develop  an  IB  method  for  the  ABLE  model.  The  description  of  IB  method  will  also  be  deferred 
to  later  because  it  is  a  large  research  and  development  problem. 

4.1  Inflow  Boundary 

The  inflow  properties  must  be  prescribed  at  the  inflow  boundary.  The  wind  components  U,  V,  W, 
and  potential  temperature  can  interpolated  either  from  the  observational  data  or  from  a  larger- 
scale  model  (e.g.,  mesoscale  numerical  prediction  model)  results.  For  an  unsteady  problem,  a 
new  inflow  velocity  value  is  prescribed  at  each  time  step.  In  the  case  of  fully  turbulent  flow,  a 
logarithmic  wind  profile  is  imposed  at  the  surface  layer.  For  example,  in  idealized,  plane-parallel 
flow  from  the  west  boundary,  the  U  wind  profile  will  be  logarithmic  while  V  and  W  are  set  to 
zero. 

Since  the  inflow  velocity  boundary  is  prescribed,  it  should  be  kept  constant  during  the  outer  or 
inner  iteration  for  a  steady  solution.  The  pressure  correction  P  is  kept  as  zero  due  to  the  constant 
inflow  velocity  in  the  steady  flow  at  the  boundary. 

4.2  Outflow  Boundary 

At  an  outflow  boundary,  the  zero  gradient  boundary  conditions  are  set  for  every  variable.  For  an 
unsteady  flow,  the  flow  variables  are  also  extrapolated  according  to  the  zero  gradient  condition 
at  the  outflow  boundary. 

4.3  Wall  Boundary 
4.3.1  Wind  and  Pressure 

For  momentum  equations,  a  non-slip  boundary  Dirichlet  condition  is  applied  at  the  walls  by 
setting  all  velocity  components  to  zero.  Besides  the  non-slip  boundary  condition,  there  is 
another  Neumann  boundary  condition  that  needs  to  be  applied  for  the  diffusive  fluxes  terms  at 
the  walls  (figure  3). 
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Figure  3.  A  schematic  diagram  for  the  wall  boundary. 


For  a  wall  boundary,  taking  a  bottom  wall  as  an  example, 


dU  dV  dW  n  dU  n  dV  n 

—  =  —  = - =  0  ,  — *0  ,  — *0. 

dx  dy  dz  dz  dz 


(32) 


The  diffusive  flux  at  the  bottom  wall  is  then  expressed  as 


™ds»#bv’-u‘ 

dz  zP  -  zB 

(33) 

dV  V  -V 

dS«pSb-Z - 

dz  zP  -  zB 

(34) 

For  this  model  discretization,  we  applied  collocated  the  computation  grid.  All  CVs  extend  to  the 
boundary  and  boundary  pressure  is  needed  to  calculate  the  pressure  forces  in  the  momentum 
equations.  A  linear  extrapolation  is  applied  to  obtain  the  pressure  at  the  boundaries. 

4.3.2  Temperature,  Thermal  Flux  from  Ground 

These  boundary  conditions  will  be  discussed  in  future  development  of  the  ABLE  model. 

4.4  Symmetry  Plane  (Free-slip) 

For  a  symmetry  plane  or  top  the  computational  domain,  a  symmetry  plane  can  be  used  to 
eliminate  momentum  transport  to  outside  of  the  domain.  In  this  type  of  boundary  condition,  the 
normal  velocity  to  the  symmetry  line  or  plane  is  zero.  Figure  4  is  an  example  of  horizontal  at  top 
of  the  boundary  layer  in  three  dimensions. 


Figure  4.  A  schematic  diagram  for  top  symmetry  boundary. 
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(35) 


dU  dV  n  dW  n 

- =  —  =  0,  - *0 

dz  dz  dz 

The  diffusive  flux  at  the  horizontal  symmetry  plane  is  computed  as 


The  boundary  conditions  are  applied  for  momentum  equations  in  such  a  way  that  the  flux 
coefficients  for  the  CV  at  the  boundary  are  added  for  the  planes  of  boundary  finite  volumes. 

4.5  Immersed  Boundary  Conditions 

Recently,  there  has  been  significant  advancement  in  IB  methods.  The  IB  method  (Tseng  and 
Ferziger,  2003;  Mittal  and  Iaccarino  2005;  Lundquist  et  al.  2007)  has  been  recognized  as  one  of 
the  most  accurate  and  efficient  methods  for  complex  geometry  in  geophysical  flows.  The  IB 
method  specifies  a  boundary  value  in  such  a  way  as  to  simulate  the  presence  of  a  body  surface  in 
a  cell  without  altering  the  structured  computation  grid  so  arbitrary  shapes  can  be  handled.  This 
method  combined  with  the  structured  box  type  grid  will  significantly  reduce  the  difficulty  in  grid 
generation  for  very  complex  lower  boundary  surface.  Because  this  is  a  topic  that  has  to  be 
studied  and  formulated  for  the  ABLE  model,  we  defer  the  detailed  description  to  the  future  once 
we  have  finished  and  tested  implementation  of  the  method. 

5.  Preliminary  Results 

It  is  a  good  practice  in  the  model  development  process  to  test  the  model  with  some  well-known 
laboratory  results  or  analytical  solution.  We  have  done  some  preliminary  tests  and  compared 
them  with  the  laboratory  results.  The  first  case  is  a  2-D  simulation  of  the  flow  past  an  infinitely 
long  cylinder,  in  which  the  well-known  lee  side  vortices  are  generated.  The  second  case  is  a  3-D 
lid-driven  cavity  flow.  Extensive  laboratory  results  are  available  from  the  literature  for  both 
cases. 

5.1  Flow  Past  a  Infinite  Long  Cylinder 

A  2-D  test  simulation  of  steady  flow  past  a  cylinder  is  performed.  The  flow  is  dependent  upon 
the  Reynolds  number,  Re  =  UtD  I  v ,  where  D  is  the  diameter  of  the  cylinder,  Ux  is  the  upstream 

fluid  velocity,  and  v  is  the  dynamic  viscosity.  There  are  many  flow  laboratory  tests  and  a 
visualization  database  is  available.  At  Re  <5  ,  the  flow  is  a  creeping  or  Stokes  flow,  in  which 
there  is  no  lee  side  vortex.  At  a  medium  Reynolds  number  around  (  R  <  50 ),  two  symmetrical 

standing  vortices  are  formed  at  the  lee  side.  At  a  higher  Reynolds  number,  these  vortices  are 
stretched  and  wavy  behavior  is  observed.  At  an  even  higher  Reynolds  number  ( Re  >  100 ),  an 
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alternating  vortex  shedding  from  the  cylinder  at  the  lee  side,  called  a  Von  Karman  vortex  street, 
is  observed. 

Preliminary  tests  for  ABLE  were  performed  for  R,,=40  and  Re=  1 00.  The  simulations  were 
performed  in  a  domain  (length  x  width  =  40D  x  20D)  large  enough  to  eliminate  the  outer 
boundary  effect  on  the  wake  formation.  The  total  grid  number  is  200  x  100.  Figure  5  shows  the 
flow  at  Re= 40,  and  R,,=  1 00,  compared  with  the  images  from  the  album  of  fluid  motion  by  Van 
Dyke  (1982).  Note  that  the  Reynolds  number  are  not  exactly  the  same  as  in  our  simulation,  but 
the  Re  are  close  enough,  i.e.,  within  the  same  regimes  of  the  flow  pattern.  The  results  also 
compared  well  with  the  numerical  simulation  of  Tseng  and  Ferziger  (2003). 

(a)  ABLE  model  Re=40  (b)  Van 


(c)  ABLE  model  Re=100  (d)  Van  Dyke  (1982)  Re=140 


Figure  5.  A  comparison  of  the  ABLE  solution  with  the  laboratory  visualization  (Adopted  from  Van  Dyke  1982)  at 
different  Reynolds  numbers. 

5.2  Three-dimensional  Lid-driven  Cavity  Flow 

This  test  case  is  specifically  for  the  3-D  framework  of  the  ABLE  model.  In  this  case,  a  cubic 
cavity  is  filled  with  a  fluid  and  the  top  cover  of  the  cavity  is  moved  at  a  steady  speed  to  specify  a 
Reynolds  number.  Several  laboratory  test  data  sets  are  available  with  different  Reynolds 
numbers.  We  use  the  test  results  by  Parasad  and  Kosseff  (1988)  for  comparison. 

In  this  simulation,  the  computational  grid  is  a  uniform  grid  with  100  x  100  x  100  grid  points.  The 
Re= 3200  was  set  for  the  comparison  with  the  laboratory  results.  For  this  simple  lid-driven  cavity 
flow,  the  flow  showed  complex  patterns.  This  result  is  the  steady-state  solution  for  the  flow.  The 
results  at  the  different  cross  sections  are  shown  in  figure  6.  The  ABFE  model  simulated  U  and  W 
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components  at  the  center  plane  at  y=l/2W  show  good  agreement  with  the  laboratory  test  results 
of  Parasad  and  Kosseff  (1989). 


ABLE  preliminary  3D  test  result  (Re=32oo) 

Lid  Driven  Cavity  flow  (3D) 


Figure  6.  ABLE  model  simulation  of  a  lid-driven  cavity  flow.  The  cross  sections  denoted  from  the  color- 
coded  planes.  The  results  were  compared  with  the  laboratory  test  data  from  Prasad  and  Kosseff 
(1989). 

We  have  also  compared  ABLE  model  results  with  other  numerical  simulations.  Figure  7  is  a  plot 
of  symmetry  plane  (Z-Y)  at  X=0.75  for  VW  components  from  simulations  of  Zang  et  al.  (1994) 
and  the  ABLE  model.  Taylor-Goertler-like  vortices  are  evident  near  the  Z=0.15.  These  types  of 
vortices  are  also  observed  in  the  laboratory  test  (Prasad  and  Kosseff,  1989).  The  flow  in  cubic 
lid-driven  cavity  is  symmetric  at  the  XY  and  YZ  planes  due  the  boundaries  and  forces  in  this 
flow.  This  property  is  also  a  good  check  for  model  code  correctness.  The  results  in  figures  6  and 
7  show  this  symmetric  property. 
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(Adopted  from  Zang  et  al.  1994) 
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Figure  7.  ABLE  simulation  compared  with  other  numerical  model  results  for  the  YZ  plane  at  X=0.75.  Left 

panel:  streamline  from  the  simulation  of  Zang  et  al.  (1994).  Right  panel:  the  flow  field  cross  section  at 
the  same  location  as  Zang  et  al. 


6.  Conclusion 


A  framework  for  a  microscale  meteorological  model,  ABLE,  has  been  developed.  A  lot  of 
difficult  tasks  are  still  ahead  of  us  to  enable  a  complete  model.  The  objective  the  ABLE  model  is 
to  simulate  microscale  (defined  meters  to  hundreds  of  meters  in  space  and  minutes  in  time) 
atmospheric  processes  in  the  ABL.  We  have  chosen  a  set  of  incompressible  set  of  governing 
equations  for  the  model  according  to  our  objectives  and  intended  scope  of  applications.  In  this 
framework,  we  have  used  a  Cartesian  coordinate  system,  a  finite  volume  approach,  and 
collocated  grids  for  the  numerical  discretization  of  the  model  system.  The  numerical  integration 
for  the  scalar  and  momentum  equations  is  described  in  detail  in  the  report.  The  SIMPLE  method 
was  selected  for  the  semi-implicit  coupling  of  the  velocity  and  pressure  fields.  Some  simple 
boundary  conditions  are  discussed  and  simple  test  cases  are  presented  in  this  report.  Some 
planned  features  for  this  model  system,  such  as  turbulence  model  for  sub-grid  scale  and  the 
treatment  of  complex  boundary,  are  also  described  briefly. 

This  is  the  starting  point  of  ABLE  model  development.  Several  planned  development 
components  of  the  system,  such  as  IB  method  for  treatment  of  complex  boundary,  turbulence 
parameterization  for  the  ABL  flow,  radiation  and  surface  energy  exchange,  and  parallelization  of 
the  computation,  will  be  carried  out  in  future  research  and  development. 
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List  of  Symbols,  Abbreviations,  and  Acronyms 


3DWF 

Three-dimensional  Wind  Field 

ADI 

alternate-direction  implicit 

ABL 

atmospheric  boundary  layer 

ABLE 

Atmospheric  Boundary  Layer  Environment 

BI-CGSTAB 

bi-conjugate  gradient  stabilized 

CV 

control  volume 

CDS 

central  difference  scheme 

FV 

finite  volume 

IB 

immersed  boundary 

LES 

large  eddy  simulation 

NWP 

numerical  weather  prediction 

RANS 

Reynolds  averaged  Navier-Stocks 

TKE 

turbulent  kinetic  energy 

SOR 

successive  over  relaxation 

SIP 

strongly  implicit  procedure 

SIMPLE 

Semi-Implicit  Method  for  Pressure-Linked  Equations 

UDS 

upstream  difference  scheme 
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